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Tumors constitute a wide family of diseases kinetically characterized by the co-presence of multiple 
spatio-temporal scales. So, tumor cells ecologically interplay with other kind of cells, e.g. en- 
dothelial cells or immune system effectors, producing and exchanging various chemical signals. As 
such, tumor growth is an ideal object of hybrid modeling where discrete stochastic processes model 
agents at low concentrations, and mean-field equations model chemical signals. In previous works 
we proposed a hybrid version of the well-known Panetta-Kirschner mean-field model of tumor cells, 
effector cells and Interleukin-2. Our hybrid model suggested -at variance of the inferences from 
its original formulation- that immune surveillance, i.e. tumor elimination by the immune system, 
may occur through a sort of side-effect of large stochastic oscillations. However, that model did 
not account that, due to both chemical transportation and cellular differentiation/division, the tumor- 
induced recruitment of immune effectors is not instantaneous but, instead, it exhibits a lag period. 
To capture this, we here integrate a mean-field equation for Interleukins-2 with a bi-dimensional de- 
layed stochastic process describing such delayed interplay. An algorithm to realize trajectories of 
the underlying stochastic process is obtained by coupling the Piecewise Deterministic Markov pro- 
cess (for the hybrid part) with a Generalized Semi-Markovian clock structure (to account for delays). 
We (i) relate tumor mass growth with delays via simulations and via parametric sensitivity analysis 
techniques, (h) we quantitatively determine probabilistic eradication times, and (Hi) we prove, in the 
oscillatory regime, the existence of a heuristic stochastic bifurcation resulting in delay-induced tumor 
eradication, which is neither predicted by the mean-field nor by the hybrid non-delayed models. 

1 Introduction 

Tumor-immune system interaction is triggered by the appearance of specific antigens - called neo- 
antigens - eventually formed by the vast number of genetic and epigenetic events characterizing tumors 
ll48l . So, the immune system may control and, in some case to eliminate, tumors [29]. This observation, 
fundamental to the so-called immune surveillance hypothesis, recently accumulated evidences |[28l . 

The competitive interaction between tumor cells and the immune system is extremely complex and, 
as such, it has multiple outcomes. So, for instance, a neoplasm may very often escape from immune 
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control, may be constrained in a oscillatory regime or, differently, a dynamic equilibrium with the tu- 
mor in a microscopic undetectable "dormant" steady-state l20l may also be established. In the oscil- 
latory regime both 'short term-small amplitude' oscillations ll39l l53l l32l l46ll and patterns of remission- 
recurrence ll50l 121 have been observed, i.e. the alternation of long dormancy phases where the immune 
surveillance is not definitive with tumor escape phases. The latter case has important and negative impli- 
cations since, on the one hand, a dormant tumor may eventually induce metastases through blood vessels 
formation and, on the other hand, the neoplasm may develop strategies to circumvent the immune system 
action, thus restarting to grow Il54ll48ll28ll5ll . This evolutionary adaptation, termed "immunoediting", 
typically happens over a significant fraction of the average host life span l28l and, among its many ef- 
fects, it negatively impacts on the effectiveness of immunotherapies ll22ll . These therapies, consisting in 
stimulating the immune system to better fight, and hopefully eradicate, a cancer, are a simple and promis- 
ing approach to the treatment of cancer ll27l . even though a huge inter-subjects variability is observed, 
which makes the results of immunotherapy clinical trials quite puzzling lffll4ll38Tl. 

As far as the modeling of tumor-immune system interplay is concerned, many mean-field models 
have appeared 0T] @3j @2j HU US HH, some of them including delays HO] [52l ESJ. However, since 
tumor cells exchange a number of chemical signals with other kind of cells, e.g endothelial cells or 
immune system effectors, they are an ideal object of hybrid modeling where some agents are represented 
by discrete stochastic processes, especially those in low numbers ll34l . and chemicals are represented by 
mean-field equations lPT2l |2~T1 . This allows to consider the intrinsic noise of the model and, when the 
mean-field approach would be an over-approximation, this may provide more informative forecasts lfl2l . 

In lfl2ll2Tl we proposed a hybrid version of the well-known Panetta-Kirschner fl4~T1 mean-field model 
of tumor cells, effector cells and Interleukins-2. The original model forecasts various kinds of experimen- 
tally observed tumor size oscillations P9ll53ll32ll46ll50l . as well as microscopic/macroscopic constant 
equilibria. However, its hybrid analogous suggests - in addiction to replicating original deterministic 
forecasts - that immune surveillance, i.e. tumor elimination by the immune system, may occur through 
a sort of side-effect of large stochastic oscillations. By discretizing both tumor and effector cellular 
populations, and by approximating the interleukins with a mean-field equation, probabilistic tumor erad- 
ication times s have been quantitatively determined for various model configurations. Also, in ||2T1 the 
model was extended to account for both interleukin-based therapies and Adoptive Cellular Immunother- 
apies, i.e. the transfusion of autologous or allogeneic T cells into tumor-bearing hosts 071 . and model 
outcomes have been investigated under various therapeutic settings . 

However, that hybrid model did not take into account that, due to both chemical transportation and 
cellular differentiation/division, the influence of tumor on immune system effectors recruitment and pro- 
liferation is not instantaneous but, instead, it exhibits a lag period. Thus, to represent this phenomenon, 
we here couple the mean-field equation for Interleukins-2 with a bi-dimensional delayed stochastic pro- 
cess describing such a delayed interplay. This delay serves to approximate missing dynamical com- 
ponents, e.g. exchanged chemical signals, maturation and activation of T-lymphocytes mediated by 
B-lymphocytes ll30l or, more in general, the fact that the immune system needs time to identify a tumor 
and react properly [49 ]. Of course, a full phenomenological model of these processes would be desirable. 
However, attempting to model each relevant stage of this process is currently impossible also because of 
the lack of systematic data iTTOl . Thus, despite this abstraction being a highly macroscopical and simplis- 
tic representation of tumor-immune system interplay, it can still provide useful insights in understanding 
this very fundamental and complex interaction. 

This new hybrid system with delay is a stochastic process combining the Piecewise Deterministic 
Markov process ll24l underlying the delay-free model (6) [71 [H with a superimposed clock structure of 
a Generalized Semi-Markov process 051 . as one of those underlying chemically reacting systems with 
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delays ifTTl [T3l . As a consequence, numerical realizations of the model are obtained by combining a 
Gillespie-like Stochastic Simulation Algorithm with delays [2] with the algorithm to simulate the delay- 
free hybrid system [12]. Via numerical analyses (i) we study the effect of various delays on tumor mass 
growth, (z7) we quantitatively determine eradication times as probability distributions, (Hi) we define a 
time-dependent sensitivity coefficient relating tumor mass and delay amplitude and (rv) we prove, in the 
oscillatory regime, the existence of a heuristic stochastic bifurcation resulting in delay-induced tumor 
eradication, which is neither predicted by the mean-field model nor by the hybrid non-delayed model. 

The paper is structured as follows. In Section |2] we present the model with delay, discuss its for- 
mulation in terms of hybrid automata and the underlying stochastic processes. In Section [3] we discuss 
algorithms for the realization of such processes and, in Section 01 we present the results of our simula- 
tions. Finally, in Section|5]we draw some conclusions and discuss future works. 



2 Model definition 

We start by extending the model given in [12, 21] with the simple form of constant delay in the immune- 
response. We consider two cell populations, i.e. tumor cells T and immune system effectors E, and 
the molecular population of Interleukins-2 (IL-2) /. A Delay Differential Equation (DDE) model can be 
stated by considering two equations for cells 

r / = i-rfi--r)- PtT e e' = -^-E-ti E E+cT{t-e) (i) 



V ) g T V + T g E +I 
and one equation for ILs-2, that is 



(2) 

V giV + T 



These equations are obtained, as in 1211 . by converting into total number of cells the densities and E* 
of the mean-field model in PTI (not shown here), i.e. T* = T/V and E* = E/V where V is the blood 
and bone marrow volumes for leukemia. In lfT2l an hybrid model is built by switching to a discrete 
representation of the populations ruled by equation (JJ) and by keeping continuous IL-2, as we shall 
discuss in the following. An immediate consequence of this is that equation £T|) is interpreted as a 
set of stochastic events, whereas equation (fSJ) is left unchanged. In this model the tumor induces the 
recruitment of the effectors at a linear rate cT(t — 8) with delay 8 > 0. With respect to l2ll . where 
instead the recruitment is instantaneous, i.e. 8=0, the delay effect is to approximate missing dynamical 
components ll30l[49l . As in the original model formulation c is a measure of the immunogenicity of the 
tumor, i.e. c is "a measure of how different the tumor is from self" BTI . Biologically, c corresponds to the 
average number of antigens, i.e. secreted antibodies and/or surface receptors on immune system T-cells, 
expressed by each tumor cell. Interleukins stimulate effectors proliferation, whose average lifespan is 
ilj; 1 , and the average degradation time for IL-2 is 1 . The source of interleukin is modeled as depending 
on both the effectors and the tumor burden. Michaelis-Menten kinetics rule IL-2 production by the tumor 
immune-system interplay, effectors recruitment by their interplay with IL-2 and effectors-induced tumour 
death. Finally, tumor growth is logistic with plateau l/b . 

In lfT2l it is shown that, when 8=0, the hybrid model predicts a desired tumor eradication via 
immune surveillance, whereas the mean-field analogous does not BTI . Subsequently, in [21] Adoptive 
Cellular Immunotherapies and Interleukin-based therapies are added to the model. By focusing on real- 
istic therapeutic settings, i.e. impulsive and piece-wise constant infusion delivery schedule, it is shown 
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that the delivery schedule deeply impacts on the therapy-induced tumor eradication time. The advantage 
of resetting the mean-field version model to the hybrid setting allows to quantitatively determine the 
probability of eradication, i.e. ZP[T(t) = 0] for some t, given various model configurations. 

In hybrid systems terminology, when 6 = this model is a Stochastic Hybrid Automaton (SHA, (6) 
7 ]) with modes in N x N, i.e. the "control" part of the automaton, recording the cellular concentrations. 
The SHA consists of a mode for each possible value of E and T, i.e. a mode q = (q E ,q T ) to count q e and 
qj effector and tumor cells, with inside the vector field of equation ((U), i.e. such a mode contains 

lit) =B q + (I q - B q ) exp {-^{t - t q )) (3) 

with initial condition l(t q ) = I q when t q is the mode entrance time and B q = [piqrqE / (giV 2 + qrV)]/lX]. 
An automata execution switches probabilistically between modes, while continuous paths of I(t) are 
determined; so, when jumping from mode q, at time t q , to mode q', at time t q i, the initial condition of 
/(?), i.e. I(t q >), is set equal to the last evaluation of I(t), i.e. I(t q ). Jumps between modes are determined 
by the time-inhomogenous stochastic events, i.e. the jump rates triggering changes in E and T depend 
on /(f) fl2j . The exit times for mode q are given by the time-dependent cumulative distribution function 

9> q \z\ = exp M£ jf a i>q (t q + t)dt\ (4) 

and the probability of jumping to mode q' , given the exit time i, is 

'LjeQ a j,q( t q + ' V 



^q[q'\z] 



if Q = {j I q + v, =q'\ 
Ziai, q (t q + r) (5) 

otherwise . 



Notice that two stochastic events, i.e. a2, q and a? A , trigger jumps to the same new mode, i.e. jumps from 
q = (qE,qr) to (q E — \,qr), so their probabilities sum up in Q. Here the Gillespie-like 11331 notation is 
used so Vj is the j'-th column of the system stoichiometry matrix 

1-1-10 
1 -1 1 

and the jump rates in mode q = (qr,qE) are the time-dependet propensity functions [34] 

a\,q{t) = r 2 qr a 2,q(t) = r 2 bV~ l q T (qT ~ 1) 

a^, q {t) = {pTqTqE)/(grV + q T ) a 4 , q (t) = \p E qEl(t)]/[gE+I{t)] 

a5,q (t) = \Ie(Ie ae, q (t) = cq T . 

Notice that all but at, q are time-homogenous jump rates, i.e. do not depend on the I(t) inside the mode, 
but, because of a^ q the underlying stochastic process is not homogenous. 

Executions of this SHA are trajectories of the underlying Piecewise Deterministic Markov Process 
(PDMP, |[24l "). a jump process over vector fields which behaves deterministically and whose jumps are 
triggered by (i) hitting user-defined boundaries of the state space and (if) time-inhoumogenous jump 
distributions. Actually, for this case, the underlying PDMP has no hitting boundaries but only time- 
dependent jump rates linked to the vector field I(t). The state space for the PDMP is N x N x R + , 
as shown in Figure [TJ In there, once the process enters state (qx,qE,qi) the only movement gradient 
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Figure 1 : State space for the hybrid model. The state space N x N x M. for the PDMP [24] underlying 
the hybrid model when 6=0. Once the process enters state (qx,qE,qi) the only movement gradient 
is on the z-axis, i.e. the horizontal component (qr,qE) is fixed and the process moves according to the 
vertical vector field represented by the empty arrows. The process persists moving according to equation 
(@]), and then moves on the N x N sub-space, i.e. the horizontal discrete grid denoted by the full arrows, 
according to equation (f5]). When 6 > the process is enriched with a clock structure as for GSMPs ||35l . 
thus inducing further jumps on the horizontal discrete grid to account for delayed reactions. 

is on the z-axis, i.e. the horizontal component (qr,qE) is fixed and the process moves according to the 
vertical vector field. The process persists moving according to equation (O, and then moves on the N x N 
sub-space, i.e. the horizontal discrete grid, according to equation ©. 

When 6 > the SHA jumps are no more given by a continuous time Markov process but, instead, 
by a Generalized Semi-Markov Process (GSMP, j35l0 . a kind of process characterizing a large class of 
discrete-event simulations ifTTl [T6l Q. It is shown in |[T3l [TTll that these process underly Gillespie- 
like lT34l chemically reacting systems with deterministic delays, those indeed used here. In these discrete 
processes (i) the embedded state process is a Markov chain and (if) the time between jumps is an arbi- 
trarily distributed random variable which may depend on the starting and the ending modes. When (a) a 
single jump event is present in each state then the process is a Semi-Markov Process, when (b) multiple 
are currently running then the process is a GSMP and, finally, when (c) the jump times are exponentially 
distributed, i.e. memoryless, then the GSMP becomes a Continuous-Time Markov Chain (CTMC). 

We recall the definition of finite-state GSMPs as in IfTTl ; the overall process will have the structure 
of the PDMP with the GSMP clock structure superimposed. We remark that, even if the state-space of 
our process is not finite, i.e. both T and E can theoretically grow unbounded, we could arbitrarily define 
two thresholds to limit the cells growth to account for biologically realistic configurations. Regions of 
the parameters in which unbounded growth of the cellular populations are determined in HT1 [121 . and 
could be used to define such thresholds. Here, since we only perform simulation-based analysis of these 
processes we can avoid restricting the GSMP to a finite state space. Let E = {e\ , . . . ,e n } be a finite set of 
events and, for any state s € S, let s t- > E (s) be a mapping from s to a non-empty subset of E denoting the 
active events in state s. In this GSMP one exponential event is always the one related to the jump process, 

1 Theoretically, this process might be equally refrained as a pure PDMP with unbounded number of clocks and infinite 
dimensional state space. Even though proving existence and uniqueness of the solutions of the ODE would be feasible, we 
think that the combined process allows for the definition of an efficient simulation algorithm (see Section[3}. 
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Algorithm 1 Input: (Tq,Eq,Iq), start time to, stop time t stop 

1: set initial mode q <— (qT ,lE ) an d set I (to) = Io', 

2: while t < t^p do 

3: let r\ ~ U[0, 1] determine the mode exit time T as ^[f] = 1/n thus solving equation ©; 
4: determine the jump rates + t), set /(? + t); 
5: jump to mode with probability HP q [q' | t]; 

6: end while 



and there is one event for each delayed transition pending; in next section an algorithm to simulate this 
joint process is given. When in state s the occurrence of one or more events triggers a state transition, 
the next state s' is chosen according to a probability distribution p(s';s,E*) where E* C E(s) is the set 
of active events which are triggering the state transition. Clocks are associated with events and, in state 
s, the clock associated with event e decays at rate r(s,e) = 1 since, in this case, time flows uniformly for 
the involved components. When, in a state s, there are no outgoing transitions, i.e. E(s) = 0, the state 
s is said to be absorbing and it models a terminating process. The set of possible clock-reading vectors 
when the state is s is 

C(s) = {c=( Cl ,...,c M ) | c i £[0,oo)Ac i >0^e i £E(s)} 

where c; is the value of the clock associated with e, ; c, G ^ where ^ is the set of clock evalutions. In 
state s with clock-reading vector c = (c\ , . . . , cm), the time to the next transition is 

t*(s,c) = min c,7r(j,e;) = min c; 

{i\e,eE( S )} 0>;&EM} 

where Cj/r(s,ei) = +°° when r(s,ej) = 0. The set of events triggering the state transition is then 

E*(s,c) = {e; EE(s) I Ci-t*(s,c)r(s,ei) = 0} . 

Actually, as is shown in lfl3l . by probabilistic arguments it is possible to show that, for chemically 
reacting systems with delays, there is a unique possible events triggering at once, i.e. E*(s,c) is a 
singleton. When a state transition from s to s' is triggered the events E* expire, leaving E'(s) =E(s)\E*. 
Moreover some new events are created; this set of new events is E(s') \ E'(s). For these events e' a clock 
value x is generated by a distribution-assignment function F(x;s' ,e' ,s,E*) such that F(0;s' ,e' ,s,E*) = 
and \im x ^ tao F(x;s l ,e' ,s,E*) = 1. For the old events in E(s') P\E'(s) the clock value in state s at the time 
when the transition was triggered is maintained in In s' events in E'(s) \ E(s') are cancelled and the 
corresponding clock value is discarded. The GSMP is a continuous-time stochastic process {X(t) \ t > 0} 
recording the state of the system as it evolves and its semantics is given in terms of a general state space 
Markov chain storing both the state of the process and the clock-reading vectors (32). 

3 Simulating the model 

We present here an algorithm to realize trajectories of the the underlying PDMP with the superimposed 
GSMP clock structure and provide model parameters. 
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Algorithm 2 Input: (Tq,Eq,Iq), start time to, stop time t stop 

1: set initial mode q <— (qT ,lE Q ), set I (to) = Io and empty scheduling list IT; 

2: while t < t^p do 

3: let r\ ~ U[0, 1] determine the mode exit time % as & q \x\ = \/r\ ; 

4: if t < head(n) then 

5: determine the rate triggering the jump according to a.j,q(t + t), set /(f + t); 

6: if the jump is triggered by ae, q then 

7: stay in mode q, set t <—t -\-T and schedule, i.e. enqueue(? + T + 6, Pi); 

8: else 

9: jump to mode q' with probability <$^[</ | t]; 

10: end if 

11: else 

12: let t' = head(n), jump to mode + l,qr), set /(? + t'), dequeue(n) and t <— t + t'; 

13: end if 

14: end while 



Model simulation. When 6=0 the SHA trajectories are generated by Algorithm [Q an extension of 
the Gillespie Stochastic Simulation Algorithm (SSA) ll33ll34l accounting for time-dependent jump rates 
and specifically tailored for this hybrid system EH . Jump times are given by solving equation (@]). 

When 6 > a combination of such an algorithm with the SSA with Delays (DSSA, (3] (HI) is 
required. The DSSA generates a statistically correct trajectory of the GSMP underlying chemically- 
reacting systems with delays ifTTl [T3l . Practically, such an algorithm is the SSA wrapped within an 
acceptance/rejection scheme to schedule/handle reactions with delays. Thus, the DSSA provides an 
algorithmic approach to the solution of the Delay Chemical Master Equation (DCME, [11]), the non- 
Markovian master equation ruling chemically reacting systems with delays. In this hybrid case, the 
system master equation is defined over the hybrid state-space [31, 14J and extended to account for the 
delays, i.e. it is a differential Chapman Kolmogorov equation with delays. 

We present here Algorithm [2] which, at the best of our knowledge, is the first attempt to combine an 
algorithm for hybrid systems with delays, in the context of biological Gillepie-like systems. This should, 
in turn, suggest further extensions towards the formal definition of SHA with delays. The algorithm uses 
a acceptance/rejection scheme and a scheduling list IT, as other DSSAs do. In this case, since a unique 
reaction with constant delay is present, IT is a standard queue data structure offering head, dequeue 
and enqueue operations. The algorithm works by determining, at each iteration, both the exit time from 
the current mode and the next mode, if any, or the scheduled reaction to handle. So, when at time t q 
the automaton enters a mode q, the exit time T (step 3) is determined by the parallel solution of I(t), 
t > t q , and & q \x\ as triggered by the jump rates aj, q (t)- As in lfl2l . samples from <^[t] are obtained by 
a unit-rate Poisson transformation (step 3), i.e. 

Y^J^ a^ q (t q + t)dt = ln 

with r\ uniformly distributed. Notice that in this equation, whose analytical solution is unknown, the 
computation is speeded up by using the analytical definition of I(t), i.e. equation ([3]). If no reactions 
with delays are scheduled to complete in [t q ,t q + t], i.e. z < peek(TT), the new mode is chosen as in the 




G.Caravagna, A.Graudenzi, A.d'Onofrio, M.Antoniotti & G.Mauri 



113 



SHA for 6 = by a weighted probabilistic choice depending on a^ q (t + t), i.e. the j satisfying 

j-l 6 ; 

£ a i;q (t + t) < r 2 £ a* i9 (f + t) < £ a /)9 (* + t) 

i=l k=l i=l 

with r2 uniformly distributed. However, if the jump is induced by the rate with delay, i.e. a^ q , the 
automata stays in mode q and the effectors recruitment is scheduled at time t q + t + 8 by means of the 
enqueue operation. This corresponds to assuming the purely delayed interpretation of delays ifTTl 121. 
being a reaction with no reactants. Finally, if a reaction with delay is scheduled in [t q ,t q + t], then the 
jump time is rejected, the system moves to the time at which the reaction is scheduled, a new effector cell 
is recruited, i.e. the system jumps from mode (qE,qr) to mode (q E + l,q T ) and the scheduled reaction 
is dequeued from IT. 

Model parameters. We use parameter values taken from lfl2ll . The baseline growth rate of the tumor 
is r = OASdays^ 1 and the organism carrying capacity is b = l/10 9 m/ _1 . The baseline strength of the 
killing rate of tumor cells by E, of the IL — 2-stimulated growth rate of E and of the production rate for 
I are, respectively, pj = Iml/days, p E = 0.l245days~ l and pi = 5 pg/days. The corresponding 50% 
reduction factors are gj = 10 5 ml , gE = 2- 10 1 pg/ 1 and gj = 10 3 mZ _1 , respectively. The degrada- 
tion rates are He = 0.03 days^ 1 for the inverse of the average lifespan of E and /I/ = 10 days~ l for the 
loss/degradation rate of IL2. Finally, the reference volume is V = 3.2ml. 

These values pertain to mice ATI l40l and are taken from ll25l l43l . where accurate fitting of real 
data concerning laboratory animals were performed. Volume V, instead, has been estimated in Ifl2l by 
considering the body weight and blood volume of a chimeric mouse. The value of 6 and c are varied in 
each configuration and given in the captions of figures. 

4 Results 

With the purpose of investigating the effect of different delays on the tumor eradication time, if any, 
and on the tumor growth size, we performed extensive simulations of various model configurations. 
All the simulations have been performed by a JAVA implementation of the model running on the cluster 
scilx . disco . unimib . it, i.e. 15 dual-core nodes, 2.0Ghz, processors and 1 GB of memory. Simulation 
times increase as T and E increase in size, spanning from few minutes to some hours, thus requiring a 
cluster capabilities to perform thousands of simulations in reasonable time. 

We always used the initial condition {Tq,Eq,Iq) = (1,0,0), one of those used in lTT2l where also 
the effect of an initial bigger tumor or effectors mass is investigated. However, we here use this initial 
condition since it allows to observe various qualitative behaviors [12]. For c = 0.02, a value used in 
Figure 2 of lfl2l . we used 6 £ {0,0.5, 1, 1.5,2,2.5,3} since, for 6 > 3, it is shown in [23] that the tumor 
mass grows up to the carrying capacity of the organism, i.e. l/b. We remark that 6 units are days, and 
6 > 3 is a biologically unrealistic value as shown in |[23l . We performed 10 3 simulations for each delay 
configuration, and we plot in Figure|2]the averages tumor and effectors growth, i.e. (T(t)) and (E(t)). 

Notice that, even though in each configuration the model still predicts tumor eradication, the tumor 
mass grows significantly more for higher delay values, i.e. for 6 = it reaches around 10 6 cells whereas 
for 6 = 3 it is 5 times bigger. This, in turn, stimulates the immune-response as shown by the plots of 
the empirical probability density of the eradication time, i.e. &[T(t) = 0] with t €f N. Notice that, even 
though the state with T = is not absorbing in the GSMP, i.e. further reactions would lead to the natural 
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Tumor cells e = n 9=1 

5 0E6 0M • ■ cc ' 




Figure 2: Tumor and Effectors growth and eradication times. In left we plot the average growth 
(T(t)) and (£(/)) as of 10 3 simulations with c = 0.02, d £ {0,0.5,1,1.5,2,2.5,3} as reported in the 
legend and (To,Eo,Io) = (1,0,0). On the ;t-axis days are represented, on the y-axis number of cells. In 
right we plot the empirical probability density of the eradication time, i.e. &[T(t) = 0] with t £ N, for 
9 £ {0, 1,2,3}. On the y-axis probability density are represented. 

death of effector cells ending to the absorbing (0,0) state, this corresponds to estimating the expected 
time for a quasi-absorbing state. These plots suggest that, though the tumor mass grows more and more 
rapidly for higher 6 - as one might expect - the effect of the consequent immune response is also larger, 
inducing a quicker eradication of the tumor, given that the mean peak for 6 = is around day 125, 
whereas for 6 £ {1,2,3} is around day 120, 118 and 115, respectively. This is a rather counterintuitive 
result, which hints at a functional role of delay in controlling the expansion of the tumor mass. 

In order to quantitatively determine the sensitivity of tumor growth with respect to 6, we perform 
parametric sensitivity analysis (PSA) by using the technique defined in |[T8l . which we now briefly 
recall. This is a numerical procedure specifically defined for discrete stochastic models; it is numerical 
since models are only rarely analytically solvable. As model output variable we use the whole ^[T(t)] 
rather than, for instance, its overall mean or mode, to capture dramatic variations in ^[T(t)], potentially 
induced by small perturbations on 6. Besides, we scan a wide range of values for 6, given that the overall 
dynamics can be differently sensitive in various regions of the parameters space. For this reason, in lfl8l 
the model sensitivity to a given parameter is defined as a function of the parameter itself. Differently 
from the mean-field case, where just T(t) could be used, the stochastic sensitivity is computed as in lfl5l 

d0> e [T(t)] 

Sd{t)= 36 (6) 

where &g[T(t)] is the probability of the tumor mass, given a value of 8. The sensitivity analysis is 
then based on a measure for discrete stochastic systems or, analogously, for the discrete part of hybrid 
systems, obeying a generic chemical master equation [36], i.e. 

d& e [T(t) = 



S T (t,0)=E[\s B (t)\}= [ 
Jn 



de 



&e[T{t)=x\dx. (7) 



The dependency of ^g[T(t)] with respect to 6 is then represented by a curve, which should be obtained 
as a function of a possibly large range of values of 6, instead of punctual perturbations. Here it is 
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Figure 3: Sensitivity analysis. In left we plot a 3-D representation of the sensitivity curves Sr(t,6), 
plotted in correspondence of each delay value of the form 6 = O.Ik with < k < 30 and k € N, for 
each t € [0,200], i.e. equation (0, as obtained by 3 x 10 5 independent simulations. In right we plot the 
sensitivity curve Sj(t) of equation (HJ. 



obtained by interpolating the points with a polynome of order D — 1, where D is the number of different 
values of delay. The model overall sensitivity coefficient, which does not depend on 6, is then 

S T (t)= I S T (t,d)dG (8) 

where the finite domain Q.6 for 6 is used. Notice that, since densities integrate to 1, the sensitivity 
coefficients do not require to be normalized as is the case for mean-field models. Also, the integral on N 
is discrete, and can be therefore represented as a summation. 

To apply this technique we performed 10 3 simulations for each delay value in Cld = {0.1 k | < k < 
30, k £ N}, thus we use 3 x 10 5 independent simulations, D = 30 and every density function is computed 
on the range [0;maxj], where maxj is the maximum observed value of T for all the values of 6, in all 
the simulations. The sensitivity function Sr,e(^,t) is then derived by integrating, for any 6, the absolute 
value of the derivative d^Pg[T(t)]/dd is evaluated in x € N and weighted by ^g[T(t) =x] according 
to equation f7]). Notice that this method does not discriminate the sign of the observed variatiorH. The 
sensitivity curves, i.e. equation (O and © are shown in Figure [3] 

One important general result is that the model sensitivity to the variation of 6 is not time-invariant, 
as shown in Figure[3j It is indeed possible to detect two intervals in which the influence is maximum, i.e. 
the intervals [10,25] and [115, 160], while in the other regions the sensitivity is essentially not relevant. 
In particular, the overall sensitivity magnitude is much larger in [115, 160], almost doubling the overall 
maximum of the first interval (right figure). This result suggests that a variation in the response time of 
the immune system can indeed influences the development of the tumor mass, but only in two specific 
conditions: (i) before that the tumor begins its expansion (i.e. first interval), either preventing or favoring 
it; (ii) after that the tumor has reached its maximum size, inducing either an enlargement or a reduction 
of the final eradication time. By looking at ST.e(d,t) (left figure) it is then possible to notice that (Hi) 
in regard to the first interval, the overall sensitivity is scarcely correlated to the specific 6, while (iv) 
the sensitivity curves corresponding to [115, 160] usually present a bell-shape, often characterized by a 



2 To perform PSA we only adopted Lagrange polynomial interpolation, even though multiple interpolation methods could 
be used and compared, e.g. spline or other non-linear interpolation techniques. 
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Figure 4: Stable oscillatory equilibria. In left we plot T{t) and E(t) for a single run with c = 0.035 
and 6 € {0, 1.5} as reported in the legend. The initial configuration is (Tq,Eq,Iq) = (1,0,0). On the 
x-axis days are represented, on the y-axis number of cells. In right we plot the phase space of the system 
restricted to T and E, and we show a stochastic switch to the null attractor for 6 = 1.5. 
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Figure 5: Heuristic stochastic bifurcation for 6 = 1.5. We plot the empirical probability density of 
the eradication time, i.e. &>[T{t) = 0] with t £ N, for c = 0.035, d = 1.5 and (T ,E ,I ) = (1,0,0) as 
evaluated by the 196 cases, out of 1000, in which the system jumps to the null attractor for T. 



unique maximum value of sensitivity, with respect to a specific 9. This suggests that a variation in can 
provoke different repercussions on the overall dynamics in distinct regions of the parameter's space. 

In order to investigate the role of delays for the system in the oscillatory regime, we performed sim- 
ulations with 0.03 < c < 0.035, a region for which both the deterministic system (i.e. Figure 2D of HTTP 
and the therapy-free hybrid model (i.e. Figure 7 of lfl21 ) predict tumor sustained/dumped oscillations. In 
Figure |4] (left) we plot the effect of delays in the oscillatory regime for c = 0.035, 6 € {0, 1.5} and initial 
configuration (7b,£rj,/o) = (1)0,0). Here we simulate the model for around 10000 days, i.e. 27 years, 
a value far beyond the life expectancy of a mouse - on which parameters are fitted - but which serves 
mainly to prove the stability of the equilibrium, if any. 

It is immediate to notice that, for 6 = 1 .5 the tumor mass does not seem to reach a small equilibrium, 
as instead it happens for the delay-free case. Indeed, in the former case the tumor mass spans between 
very low values and 3 x 10 5 , in the latter the oscillations are dumped up to around 10 5 cells. Furthermore, 
the first oscillation peak is around 4.5 x 10 5 for 6 = 1.5 which is a considerably bigger values than that 
one reached for 8 = 0. These amplified oscillations often arise when models are enriched with delays 
B7ll45l l441 and reach very small values as shown in Figure [4] (right) where the phase space of the system 
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Figure 6: Mean-field model. We plot deterministic simulations of model (Q]|2]) for c = 0.035, 6 G 
{0,0.5,1,1.5,2,2.5,3} (higher peaks for higher values of 0) and (r ,£b,/ ) = (1,0,0), T(t) = for 
t < 0. Notice the tumor resting period in t G [120, 160] (right zoom for 6 G {1, 1.5,2}), the length of 
which depends on 6, is the one in witch the hybrid system probabilistically switches to the null attractor 
for T. On the .x-axis days are represented, on the y-axis number of cells. 

restricted to T and E is represented and a stochastic switch to the null attractor for 6 = 1.5 is shown. 
Surprisingly, this result in some simulations showing eradication for 6 = 1.5, an unexpected outcome 
for the oscillatory regime since for d = none of 1000 simulations have shown eradication (not shown 
here). Instead, 196 out of 1000 simulations, i.e. almost 20% of the cases, for 8 = 1.5 show eradication 
reached immediately after the first spike of the oscillations. This clearly suggests the existence of a 
heuristic stochastic bifurcation close to 6 = 1.5 with a switch to the null attractor for T, i.e. T — >• 0, so 
that, for some cases, the tumor gets eradicated. In Figure [5] we plot the empirical probability density of 
the eradication time, i.e. ^[T(t) = 0] with t G N, as evaluated by these 196 cases. This conclusion is 
strengthened by observing that, for 8 G {2, 3}, the tumor is always eradicated (in 1000 cases, not shown). 

Moreover, this is an interesting outcome as compared against the predictions of the mean-field model. 
In fact, in Figure [6] we show deterministic simulations of model CD© for 6 G {0,0.5, 1, 1.5,2,2.5,3}, 
restricted to t G [0,400] and with extended analogous initial condition 

In there it is possible to observe a tumor resting period for t G [120, 160], the length of which depends 
on 6. Small values in such period are predicted, i.e. for 6 = 2 we observe T(t) < 1 and for 6 = 1.5 
we observe T(t) sa 10 in accordance with the simulations we performed. In this same period, instead, 
the hybrid system probabilistically switches to the null attractor for T, thus suggesting the importance of 
resetting the model in the hybrid setting which, as in Ifl2l . is again proved to be more informative. 

5 Conclusions 

In this paper we study the effect of a constant time delay in effectors recruitment in a tumor-immune 
system interplay hybrid model. The model, analogous of a well-known mean-field model [41], was 
proved to be more informative to forecast onco-suppression by the immune system lfl2l as a conjunction 
of the intrinsic tendency of the immune system to oscillate, significantly evidenced by the deterministic 
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model, with the intrinsic noise. This phenomen, which is triggered by the appearance of specific neo- 
antigens resulting from genetic and epigenetic events characterizing tumor cells [48 ], is fundamental to 
the immune surveillance hypothesis, a promising approach to the treatment of cancer EBl . 

Modeling such an interplay requires considering biological entities at multiple scales. As such, tumor 
growth is an ideal object of hybrid modeling 021. Extending the model in lPT2l with delays allows to 
account that, due to both chemical transportation and cellular differentiation/division, the influence of 
tumor on effectors recruitment and proliferation exhibits a lag period. Of course, an explicit model of 
the missing dynamical components, e.g. chemical signals, maturation and activation of T-lymphocytes, 
would be desirable but is currently unfeasible, also because of the lack of systematic data [10]. 

In this paper we contextualized this model within Stochastic Hybrid Automata, when the delay is 
0, so to give it a semantics in terms of Piecewise Deterministic Markov Processes ll24l . When delays 
are present we combine the underlying process with a clock structure for a Generalized Semi-Markov 
process IT351 . as for chemically reacting systems with delays ifTTTl . We present a novel algorithm to 
simulate this extended hybrid model and, via numerical analyses, we quantitatively determined the effects 
of various delays on tumor mass growth and determine the eradication times as probability distributions, 
under various configurations. Under these configurations we adopted a parametric sensitivity analysis 
technique to relate the tumor growth to the delay amplitude. Also, we have shown that the stochastic 
effects driving the system to the eradication can unexpectedly appear even in the oscillatory regime. In 
fact, in there we proved the existence of a heuristic stochastic bifurcation, which is neither predicted 
by the mean-field model nor by the hybrid non-delayed model. Thus, despite our model being a highly 
macroscopical and simplistic representation of the tumor-immune system interplay, we have shown that 
it can provide useful insights on the multitude of possible outcomes of this very fundamental and complex 
interaction, e.g. neoplasm evasion from immune control, immune surveillance and (dumped) oscillations. 

As far as future works are concerned, a further combination of this model with the immunotherapies 
studied in [21 j would be interesting. Also, the model itself could be extended so, for instance, the linear 
antigenic effect cT(t — 8) due to the tumor size could be corrected by assuming a delayed saturating 
stimulation. Similarly, the assumption that E' linearly depends on E could be corrected, as there are 
cases where this dependence might be non-linear, as outlined in |fl9l . Moreover, more complex form 
of delays could be considered, along the line of those used in mean-field models 11231 , e.g. weak/strong 
kernels. Finally, the mathematical formalization of hybrid automata with delays seems missing, thus 
suggesting possible extensions to the hybrid automata theory, along with their analysis techniques. 
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